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Abstract 

The pressure components of "soft" disks in a two dimensional narrow channel are analyzed in the 
dilute gas regime using the Mayer cluster expansion and molecular dynamics. Channels with either 
periodic or reflecting boundaries are considered. It is found that when the two-body potential, u(r), 
is singular at some distance ro, the dependence of the pressure components on the channel width 
exhibits a singularity at one or more channel widths which are simply related to ro. In channels 
with periodic boundary conditions and for potentials which are discontinuous at ro, the transverse 
and longitudinal pressure components exhibit a 1/2 and 3/2 singularity, respectively. Continuous 
potentials with a power law singularity result in weaker singularities of the pressure components. 
In channels with reflecting boundary conditions the singularities are found to be weaker than those 
corresponding to periodic boundaries. 
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I. INTRODUCTION 



The thermodynamic and dynamical properties of particles in restricted geometries are of 
great interest. They have been extensively studied in the context ofporous media l[ [2], [], Q], 
transport through narrow channels such as carbon nanotubes j^, [sj and pores in biological 
membranes [?| as well as in numerous other systems Qj. Perhaps the simplest and most 
convenient theoretical approach for studying fluids in cavities models the fluid by hard 



spheres. 
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his approach has been applied in a large number of studies (see for example, 
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131 ] ) . Recent studies of a dilute gas of hard disks in a narrow two dimensional 



channel have shown that the system exhibits a singularity in the pressure at a channel 
width equal to twice the diameter of the disks Q, Q. This is a consequence of the fact 
that the volume of the phase space available to the disks exhibits a singularity at this 
width. In particular, it has been found that for a channel with periodic boundary condition, 
the transverse component of the pressure exhibits a 3/2 singularity, while the longitudinal 
component exhibits a 5/2 singularity 14j- In the case of a channel with reflecting boundaries, 
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weaker singularities for both pressure components were found 

In this paper we extend these studies to consider a gas of "soft" disks in a narrow channel 
at low density. We consider several classes of two body potentials, with both periodic and 
reflecting channel boundaries. Our analysis shows that the pressure components are singu- 
lar at some channel widths whenever the interaction potential between two disks, u(r), is 
singular at some distance r . In particular, in the case of periodic boundary conditions and 
for potentials which are discontinuous at some r , the dependence of the transverse compo- 
nent of the pressure on the channel width exhibits a 1/2 singularity, while the longitudinal 
component exhibits a 3/2 singularity at some channel widths, which are simply related to 
r . The nature of these singularities becomes weaker for interaction potentials which are 
continuous, but they still display a power law singularity at r . The singularities of the 
case of reflecting boundary conditions are found to be weaker than those corresponding to 
periodic boundaries. Although we have not analyzed narrow channels in three dimensions, 
we expect similar phenomena to take place there as well. The nature of the singularities in 
the pressure in three dimensions is expected to be different from that of the two dimensional 
case. 

In the following sections we study several classes of two body potentials u(r) for both 
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periodic and reflecting boundary conditions. In Section [TT] we present the general formulation 
of the tools used in this study, the Mayer cluster expansion for gases at low densities, and 
molecular dynamics. In Section UTTl we analyze the case of a channel with periodic boundary 
conditions, for one-step and two-steps potentials. We also study a smooth potential which 
vanishes with a power law at a critical distance. In Section [IV] we study a channel with 
reflecting boundaries for the cases of soft disks and soft disks with a hard core. Finally, a 
brief summary is given in Section [V] 



II. GENERAL FORMULATION 

We consider N disks of diameter d and mass m, interacting via a two body potential u(r) 
at temperature T. The disks are restricted to move in a channel of length L x and width L y 
with L y <C L x . In this study we analyze the pressure components of this gas using virial 
expansion to second order in the density. We also carry out molecular dynamics simulations 
of this system. In the present study the channel width, L y , is taken to be finite throughout 
the calculation. Therefore the free energy is not extensive in L y and, thus, Euler's relation 
does not hold, namely, — PV ^ E — TS — fiN. Thus, the pressure has to be calculated by 
taking the appropriate derivative of the free energy. 

In the grand canonical ensemble the free energy of a system of N disks is given by 

F(T,V,N) = -kT\nC(T,V,z) + kTN\nz , 

N = z^lnC(T,V,z) . (1) 

Here £ is the grand partition sum, z is the fugacity, k is the Boltzmann constant and 
V = L x L y . The pressure components are evaluated by taking the appropriate derivatives of 
the free energy 

dF 



Pxx V L x 



dL x 
dF 



P V — — T 

P = \{P*x + Pyy) ■ (2) 
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To second order in the fugacity z, the Mayer expansion yields 

ln£ = ^(b lZ + b 2 z 2 ) 
1 N 1 „ , 

" = 7 = p(M + 26 2 z 2 ), (3) 

where A = h/ \j2i\mkT is the average thermal wavelength and h is Planck's constant. The 
coefficients of the expansion satisfy 

h = 1 

h = ^ 2 q{L y ) , (4) 

with 

q(L y ) = J f 12 d 2 r 12 . (5) 

Here, fi 2 = e _/3u ^ ri2 ^ — 1 is the Mayer function, and j3 = 1/kT. Using Eqs. (Ti][5]) we find 
that to order 1/v the free energy is given by 

_ 1 _lihl_i nv + 2 \nX , (6) 



kTN 2v 
from which the components of the pressure tensor are obtained: 



PxxV j _ q(L y ) ^ 



kT 2v 

\ J, Anil, \ 

(8) 



PyyV _ _ q(Ly) Ly ^(L y 



kT 2v 2v dL y 

All theoretical considerations are augmented by computer simulations. For our two- 
dimensional systems, the temperature T is computed from 

(K) = (N - |) kT, (9) 

where K is the kinetic energy and the bracket denotes a time average. Here, g is the number 
of macroscopic conservation laws, which differs for the periodic boundaries of Section IHH 
(g = 3; the total momentum is constant and is taken to vanish), and for the reflecting 
boundaries used in Section [TV] (g = 1). The diagonal elements of the pressure tensor, P aa 
for a G {x,y}, are evaluated from the virial theorem, 

P aa V = (K) + W aa . (10) 
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For impulsive interactions [16], the potential contribution, W aa , is given by 

^« = ^ EISA'S. en) 



T 

c 



Here, the sum is over all collisional events c, which instantaneously change the potential 
energy during the averaging time r, = — r£j is the a-component of the separation 

(c) 

vector of the two particles involved in the event, and &v a i denotes the velocity change for 
particle i parallel to a due to that event (the velocity change for particle j being just the 
opposite). For the continuous potentials of Section IHI CI the potential contribution becomes 

1 f T . 

W aa = - dt^2^2r a)ij f atij , (12) 

T ^° i j>i 

where f a ^j is the a-component of the force exerted on i by particle j. In all figures compar- 
ing experimental with theoretical pressures, the experimental dots represent W aa / (K). The 
theoretical smooth curves represent (P aa v/kT) — 1, where the temperature required for this 
computation is taken from Eq. (jHJ). As usual, temperature units are used, for which Boltz- 
mann's constant k is unity. Because of the equivalence of the canonical and microcanonical 
ensembles, the experimental and theoretical pressures should agree up to a term of order 
0(1/ N). To make this correction insignificant, at least 60 particles, or even 120 in many 
cases, were used for the simulations. 



An event-driven algorithm is used [16l . 1 1 71 ] for the discontinuous potentials with periodic 
(Sec. IIIII) or reflecting (Sec. fTVT) boundary conditions, for which instantaneous potential 
energy changes and boundary crossings of a particle are considered as events. For the 
continuous power-law potentials of Sec. IIII CI a hybrid code is used, which will be described 
there in more detail. 

III. NARROW CHANNELS WITH PERIODIC BOUNDARY CONDITIONS 
A. Positive step potential 

We proceed by considering the pressure in the case of a two-body step potential 

. u r < d 

u(r) = { ~ (13) 

r > d, 



where u > is a constant (in our previous paper [14J we considered the hard-disk case 
u = oo). The case u < is pathological, since the particles may collapse to form a cluster, 
as long as there is no repulsive interaction at short distances. The case of a two step potential 
with repulsion at short distances and attraction at larger distances will be considered in the 
following sub-section. Here we limit ourselves to the repulsive potential case u > 0. To 
evaluate the pressure we associate with each particle an interaction disk of radius d centered 
at its position. Two particles i and j interact with each other, if the center of j is within 
the interaction disk of i, and vice versa. In the case L y > 2d the cross sections of a particle 
(a disk of radius d) and that of its image resulting from the periodic boundary conditions 
in the y direction do not overlap. Hence the integral (jSJ) simply yields 

q(L y ) = 7rd 2 {e-' 3u - 1) for L y > 2d . (14) 

For L y < 2d we note that the area of overlap between the interaction disk of a particle and 
that of its image translated in the y direction is given by 

5(7?) = d 2 (n -2d- sin 2d) , (15) 

where i? satisfies (see Fig. flTJ)) 

L y = 2dsmd. (16) 

In this case the integral (JSJ) yields 

q(L y ) = (nd 2 - 2S(e~P u - 1) + S(e~ 2f3u - 1) . (17) 

Using this result for q(L y ), and noting that 

dS 



dL n , 



-\/4d 2 -Ll , (18) 



it is straightforward to derive the expressions for the pressure. We find that to order 1/v 
and for L y > 2d one has P xx = P yy = P with 

7^ = 1-^, (19) 

kT 2v V ' 

where q is independent of L y and is given by (JHj). On the other hand, for L y < 2d one finds 

kT ~ 1 2v ' 1 Uj 

P ~w = 1 ~^-fev /i ^( 1 - e ^) 2 • (21) 

S = l-^-^A/^^a-e-^) 2 , (22) 

fcT 2v 4w V y v > ' v 7 



FIG. 1: The interaction disks arrangement for L y < 2d with periodic boundary conditions. The 
overlap area of the two disks is S. 

where q(L y ) is given by ffTTj) . It is evident that P yy exhibits a square-root singularity at L y = 
2d as in the case of hard disks 14j. This singularity originates from the term dq(L y )/dL y 



in (jHJ). On the other hand, P xx exhibits a weaker singularity with a singular term which 
vanishes as (2d — L y ) 3 ^ 2 . The reason is that unlike the P yy component, here the singularity 
originates from q(L y ) and not from its derivative. Clearly, the pressure P, which is the 
average of the two components, exhibits a square-root singularity as the more singular P yy 
component. 

In Fig. [2] the theoretical expressions for the singularity at L y = 2d = 2 are compared 
to computer simulation results for various potential step sizes u as indicated by the labels. 
Reduced units are used for which the particle diameter d and the total energy per particle, 
E/N, are unity. The energy E/N is almost exclusively kinetic in nature with a time- averaged 
kinetic temperature T = 0.986 for u = 0.5 (top), T = 0.994 for u = 1.0 (middle), and 
T = 0.996 for u = 2 (bottom figure). These temperatures vary slightly, but insignificantly, 
with the channel width L y . The density is kept constant, N/V = 0.01. As for the case 
of hard spheres (u = oo) treated already in Ref. [ijj], the agreement between theory and 
simulation results is very satisfactory. 

As the width of the channel further decreases we expect to have more singularities to 
take place at L y = d, 2d/3, d/2, . . . , which result from the overlap of an interaction disk 
with those of further neighboring disks in the y direction. Let us analyze, for example, the 
second singularity in the pressure curve, which takes place at L y = d. For d/2 < L y < d 
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FIG. 2: Comparison of the channel-width dependence of the theoretical pressures (lines) with 
numerical simulations (points) for a system of 60 disks interacting with the step potential of Eq. 
(|13p . The density is kept constant, N/V = 0.01. The potential varies from u = 0.5 (top) to u = 2 



(bottom). The singularities appear at L y = 2 and L y = 1 as explained in the main text. Reduced 

units are used, for which d and the total energy per particle, E/N, are unity. 
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FIG. 3: The interaction disks arrangement for L y < d with periodic boundary conditions. The 
overlap area of the disk and its next nearest neighbor in the y direction is Sq. 

the disk of a particle has some overlap with its nearest and next nearest neighbors images 
which result from the periodic boundary conditions in the y direction. In order to evaluate 
the overlap integral (J5j) we note that (see Fig. [3]) 

L y = (isin^o • (23) 

The overlap area between an interaction disk of a particle and that of its next nearest 
neighbor in the y direction, So, is given by 

S = d 2 (7r- 2^o - sin 2tf ) • (24) 

The overlap integral is thus expressed as 

q(L y ) = (7rd 2 -2S + S ){e-P u - 1) + (S - 2S )(e~ 2f3u - 1) + S (e- 3 ^ - 1) . (25) 

Using this expression for the overlap integral, the pressure components can readily be cal- 
culated to yield 

PxxV _ _ q(L y ) 

kT ~ 1 2v [Zb) 

P = \{P XX + Pyy) ■ (28) 

As is also shown in Figure [2J these expressions for the singularity at L y — d — 1 compare 
very well with simulation results. As before, reduced units are used for which d and E/N are 



unity. Note that as long as u > the coefficient of the singular term ^Jd? — L 2 y is positive, 
resulting in a positive compressibility just below L y = d. 



B. Two-step potential 



In order to analyze the case of disks with an attractive potential, one has to add a repulsive 
interaction at short distances to prevent the collapse of the system into a macroscopic cluster. 
We thus consider in this section a two-step potential 

r < d , 

u(r) = { u 2 d<r <D , (29) 
r > D . 

where u\ > represents a repulsive interaction, w 2 could be either positive or negative, and 
D > d is the outer radius of u%- To evaluate the pressure we associate with each particle two 
concentric interaction disks, one with radius d and the other with radius D. Two particles 
only interact with each other, if the center of the second particle lies within the interaction 
disks of the first. It is easy to see that the degree of overlap of the disks of a particle and 
those of its nearest neighbor image resulting from the periodic boundary condition in the y 
direction are singular at L y = 2D, d + D, and 2d. Thus, the pressure curve is expected to 
be singular at these three values of the channel's width. 

We now analyze the pressure curve in more detail and consider first the upper singularity 
at L y = 2D. For L y > 2D the disks of a particle and those of its image do not overlap. 
Thus the integral fl5J) yields 

q(L y ) = n(D 2 - d 2 ) (e^ 2 - 1) + nd 2 ( e - (Sui - 1) for L y > 2D , (30) 

and q is independent of L y . As in the case of a single step potential, one finds that to leading 
order in 1/v the pressure tensor is isotropic, P xx = P yy = P, with 

P^ = l-±. (31) 
kT 2v K ' 

For d + D < L y < 2D, however, the outer disks of a particle and its periodic image overlap. 
As in the case of the single step potential, the overlap area S is given by 

S = D 2 (n -2d -sin2tf) , (32) 

where $ satisfies (see Fig. ([!])) 

L y = 2Ds\n$. (33) 
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FIG. 4: The interaction disks arrangement for the case of a two-step potential and for 2d < L y < 
d + D with periodic boundary conditions. The overlap area of the outer disks of two nearest 
neighbors in the y direction is S\, while the overlap area of the outer disk with the inner one of its 
nearest neighbor in the L y direction is Si. The angles $1 and $2 are indicated . 

The resulting overlap integral (|3j) for d + D < L y < 2D is given by 

q(L v ) = (nD 2 - nd 2 - 2S)(e-? U2 - 1) + ^(e"^ 1 - 1) + S{e~ 2Pu * - 1) . (34) 
The pressure tensor in this regime is thus found to be 



PxxV 1 _ q{L y ) 
kT ~ 2v 
P yy v q(L y ) L % 



(35) 



1 - " ttV 4 ^ 2 - m - e '" U2 ) 2 , ( 36 ) 



kT ' 2v 2v V ' yy 

where q(L y ) is given by (1541 . As in the case of a single step potential, P yy exhibits a square 
root singularity, while P xx behaves more smoothly, with a weaker 3/2 power behavior at the 
transition. The compressibility below the transition is negative. 

Finally, we consider the regime 2d < L y < d + D. In this case the outer disk of a particle 
partially overlaps not only with the outer disk of its periodic image but also with the inner 
one (see Fig. (j4])). The overlap area Si between the outer and the inner disks is given by 

S 1 = id 2 (vr-2tfi) -rf 2 sin^icos^i + ]^D 2 {tt - 2d 2 ) - D 2 sin$ 2 costf 2 (38) 
where $1 and i9 2 satisfy (see Fig. (|3J)) 

L y = d sin $1 + D sin $2 , and d cos d\ = D cos $2 , (39) 
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and, hence, 



sin ^ = 2dT y (L y ~° 2 + rf2) ' sin ^ 2 = 2DT v {L y + L)2 ~ d2) ■ (40) 

It is straightforward to express the overlap integral (jSJ) in terms of the overlap areas S and 
5*i as 

q{L y ) = (nD 2 - nd 2 - 2S + 25i) (e^" 2 - 1) + (nd 2 - 2S 1 )(e~ f3ui - 1) 

+ (S - 2S 1 ){e~ 2/3u2 - 1) + 2S 1 {e-^ Ul+u ^ - 1) for 2d < L y < d + D . (41) 



According to Eqs. ([7j) and (jSj), the pressure components may be expressed as 



kT ~ 2v 



p yy v = i _ Q(W) + - e -/3u 2 y_^_ + ht{^-^ _ e - /3 " 1 )(l - c ~ 0U2 ) dSl 
kT 2v 2v dL y v dL y 

where 



(42) 
(43) 



A Q 

— - = -2d cos(A) = -2D cos(tf 2 ) . (44) 

Vbi-Jy 

Using Eq. (1391) and (HOI) , the channel- width dependence of the pressure tensor components 
for 2d < L y < d + D is finally obtained, 

kT 2v ' y ' 

5^ = 1 _ _ ^^ 4 £,2_ L 2 (1 _ e -^2)2 _ I( e -/3n 2 _ e -/3ui) (1 _ e-ZJua) 

x(4Lj£> 2 + AL 2 y d 2 + AD 2 d 2 - L y - D A - d 4 ) 1 / 2 , (46) 

where <?(-£/,,) is given by (l4"Ij) . 

In order to establish the nature of the singularity at the transition point d + D, we 
expand dS\jdL y in Eq. (j44H in terms of the small dimensionless offset e = (d + D — L y )/d. 
Introducing the small angles 

0i = | - Mi , 02 = | - M 2 , (47) 
which, according to Eq. ( l39l) . are related to e by 

^ 2 =(^) £ ' <^ 2 = TO £ ' (48) 

we finally obtain 



dSi , / 8D n 

—L = -d x - Jl 49 

dL y \d + D v K ' 
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FIG. 5: Channel-width dependence of the pressures for the two-step potential case with the fol- 
lowing set of parameters: d = 1, D = 1.5, N = 60, N/V = 0.01 and E/N = 1. The top panel 
corresponds to a potential with two positive steps, u\ = 2 and u-i = 1, whereas the lower panel is 
for an attractive outer shell with u\ = 2 and U2 = —1. 

It is readily seen that dS\jdL y and, hence, P yy , exhibit a square root singularity. As in 
the case of a single step potential, the singularity in P xx originates from S\ rather than 
from its derivative. Hence P xx exhibits a weaker 3/2 singularity at L y — D + d. It is 
interesting to note that, depending on the values of the interaction parameters U\ and u 2 , 
the compressibility just below the transition could be either positive or negative. 

In Fig. [S]we compare the respective theoretical expressions - Eq. (I3TI) for L y > 2D, Eqs. 
(JBE1 - [37]) for d + D < L y < 2D, and Eqs. gEl USD f o r 2d < L y < d + D - to computer 
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v 

FIG. 6: Potentials for A = 10 (smooth curves) and A = 100 (dashed curves) for various values of 
K" MS indicated by the labels. 

simulation results (dots) for N = 60 particles in a narrow channel of width L y with periodic 
boundaries both in x and y directions. We use reduced units for which E/N and d are unity. 
The outer diameter D = 1.5. In the top panel, results for a two-step potential with u% — 2 
and U2 — 1 are shown. The lower panel corresponds to a true square well potential with 
u\ = 2 and U2 = — 1. The agreement is very good in all cases. 

A similar analysis may be carried out near the third singularity which takes place at 
L v = 2d. 



C. Power-law potential 

We now consider a soft potential which vanishes (continuously) at the disk boundary, 

(a(i-(^) 2 Y r<d 
u(r) = { V Kd) > ~ (50) 

[ r> d 

where the two parameters A > and k > are constants. In Fig. [6] we show a few of 
such potentials for A = 10E (smooth lines) and A = 100E Q (dashed lines) for various k 
as indicated by the labels. Here E = E/N = Y^iPi/^ m + J2i J2j>i u ( r ij) /N is the total 
energy per particle. For our numerical work we use reduced units, for which the diameter 
d, the particle mass m, and E are unity. 
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Let us analyze the nature of the singularity of the pressure components slightly below 
2d. The singularity arises from integrating the Mayer function / in the overlap region of 
the two disks. For e = (2d — L y )/d <C 1, the function / is small in this region and may thus 
be expanded in powers of u(r). To second order in u, f(r) ~ —j3u(r) + (l/2)(/3u(r)) 2 . The 
singularity in the integral §5§ arises from the non-linear term in /, which is of the order e 2K 
in the overlap region. Since according to (JTSJ) the area of this region scales as e 3 / 2 for small 
e, the singular contribution to the integral (j3j), and hence to P xx , scales as e 2K + 3 / 2 . On the 
other hand, the pressure P and its P yy component scale as e 2K+1 / 2 . Thus, for small enough 
e we expect 

P = Cl [l-c 2 (2d-L y ) 2 ^ 2 ] , (51) 

and similarly for P yy , where C\,C2 are constants. In the scaling form for P xx , the exponent 
is 2k + 3/2, and the singularity is weaker. 

To test this scaling form, we carried out numerical simulations of the model. In selecting 
the parameters A and k most appropriate for numerical simulations, one should take into 
account two competing trends. On the one hand, the singular part of the pressure is expected 
to be more pronounced for large amplitude A and small exponent k. On the other hand, as 
we argue below, the channel-width interval, where the scaling form ( 15TI) is expected to hold, 
is larger for small A and large k. Thus, in order to observe the scaling behavior one has to 
choose intermediate values of these two parameters. 

To estimate the scaling interval L y%m m < L y < 2d over which the scaling form (l5~Tj) is 
expected to hold, we note that during a typical collision two particles penetrate each other 
up to a depth 5 = d — r , where r is estimated from u(r ) = E Q , r = d\Jl — (Eq/A) 1 I k . 
The expansion of the Mayer function to second order in u fails, if the third-order term 
starts to contribute more than, say, 10 %. This failure only happens for particle separations 
smaller than r , the typical separation at maximum penetration, and, hence, for untypical 
high-energetic collisions. For typical energies and penetrations, particles will pass each other 
in the channel and contribute to the pressure scaling, if the thermally possible penetration 
depth 5 exceeds the interaction-disk overlap de due to the periodic boundaries. The upper 
bound for e is thus estimated to be e max = 5/d = (d — r ) / d, and the minimum channel 
width for which scaling is expected to hold becomes L Vtmin = 2d — de max = d + r . Thus, 
the scaling interval decreases with A and increases with k. We find that k = 2 and A = 10 
are a suitable choice, which gives a reasonable scaling range, and we consider this case first. 
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scaling with the power 4.5. For e > 0.07 corresponding to L y < 1.86a, the scaling breaks 



Note that the choice k > 1 also offers the slight numerical advantage that the particle force 
is continuous and vanishes at r = d. 

The simulation results for the pressures with potential parameters A = 10E and k = 2 
are shown by the dots in the bottom panel of Fig. [7J In the simulation we used 20 particles 
at a density N/V = 0.01. The estimated scaling interval, L ymin « 1.83d < L y < 2d, is 
indicated by the shaded area. The smooth lines are a fit of Eq. (|5T|) to the numerical data 
points for P and P yy in that range. It shows that our estimate is rather conservative, since 
the fits represent the data points reasonably well in a slightly wider interval 1.76d < L y < 2d. 
The kinetic energy per particle is about 0.998i?o and varies only marginally with L y . The 
scaling is more convincingly demonstrated in Fig. [HI where the e-dependence of the singular 

for P yy is shown. The straight line indicates the expected 
7 corresponding to L y 
down as expected. 

Next, we consider a potential with A = lOO-Eo, k — 2, which is shown in Fig. [6], and 
which resembles more realistic repulsive potentials. The channel-width dependence of the 
pressures is shown in the middle panel of Fig. UJ The estimated scaling range is much 
narrower than before, L y min = 1.95d, and is indicated by the shaded area. The smooth 
lines are fits of Eq. (j5Tl) to the data points for P and P yy in that interval and confirm that 
our scaling-range estimate is rather conservative. The kinetic energy per particle is around 
0.999£'o in this case, and varies marginally with L y . 

Finally, we consider the limiting case of a rather steep potential such as for A = 400-E'o 
and k = 2 (much steeper than the k = 2 curve in Fig. [6]), which already resembles that of 
hard disks and, therefore, should give a pressure variation with L y similar to that found in 
Ref. [14 1 . The results are shown in the top panel of Fig. The average kinetic energy per 
particle is 0.9997-E'o- For L y < 1.95<i the pressure curves are indeed very similar to those of 
a hard-disk gas of 20 disks at the same density and at unit kinetic energy per particle as is 
shown in Fig. 3 of Ref. Jjj]. Differences appear for channel widths very close to 2d, which 
are due to the expected scaling. The estimated scaling range is very narrow, L^mm = 1.975d 
as indicated by the shaded area in the top panel of Fig. But a fit of Eq. (I5TT) in the range 
1.96d < L y < 2d represents the data points for P and P yy reasonably well in that range as 
is shown by the smooth lines. 

Before closing this section, we provide some details about the molecular dynamics simu- 
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FIG. 7: Simulation results for the pressures as a function of the channel width L y with periodic 
boundaries. The potential parameter A varies from A = 400 (top) to A = 100 (middle) and A = 10 
(bottom), and k = 2. In this figure N/V = 0.01, and iV = 20. The shaded areas indicate (very 
conservative) estimates of the scaling regimes. The smooth lines are a fit of Eq. (151j) for P and 
P yy to the data in the shaded regime. Reduced units are used for which d and Eq are unity. 
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FIG. 8: Scaling of the transverse pressure P yy for the power-law potential model with A = 10 and 
k = 2 below the critical channel width L y = 2. Here, e = (2 — L y ). The slope of the straight 
line corresponds to the theoretically expected scaling, k + (1/2) = 4.5. Reduced units are used as 
explained in the main text. 

lations. They were carried out with a hybrid code combining the advantages of the event- 
driven algorithm for hard particles during the forceless streaming stage with the simplicity 
of a time-stepping integration scheme during the collision of two or more particles. The 
beginning of each pair collision was determined as in the event-driven algorithms of the pre- 
vious sections. During the collisions the equations of motion were integrated with a fourth 
order Runge-Kutta scheme. The end of each pair collision was determined by interpolation 
with a spatial uncertainty of less than 10~ 8 reduced units. The moment the last interacting 
particles separate, another streaming move is initiated. This method is particularly suited 
for low densities. It even allows to accurately follow the trajectory for models with discon- 
tinuous forces. Periodic boundaries are used. In most cases a trajectory was followed for 
two million reduced time units ^md 2 /E . 
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IV. CHANNELS WITH REFLECTING BOUNDARY CONDITIONS 



A. Soft disks: single step potential 

In this section we calculate the pressure components of soft disks in a narrow rectangular 
box with elastic reflecting boundary conditions in the y direction. Since we are interested 
in the narrow channel limit where the length of the box is much larger than the width, the 
system is not sensitive to the boundary conditions in the x direction. For simplicity we 
take periodic boundary condition in this direction. We consider disks of diameter d, which 
interact with each other via the square well potential of Eq. (1131) . but with an additional 
^-function at the center of the particles: 

(5(r) + u r < d , 
{> (52) 
r > d , 

This 5 function does not affect the particle-particle interactions, but it is responsible for the 
elastic reflections from the boundary, which confine the disk centers to the volume V = L x L y , 
where L y is referred to as the channel width. 

As a result of the reflecting boundary conditions, the system is no longer translation- 
ally invariant, and the overlap integral (j^D corresponding to the second virial coefficient is 
replaced by 

q(L y ) = j-J f 12 d 2 n d 2 r 2 . (53) 

As was done in the case of periodic boundary conditions, with each particle we associate an 
interaction disk with a radius d. Two particles interact with each other only if the center of 
a particle is within the interaction disk of the other. 

In the case L y > d the overlap integral is given by (see the top panel of Fig. [S]) 



Q( L v) = T~ 

L 



7rd 2 L y -2 [ S($ } d)dy 
Jo 



d 

fill 



e-l» - 1) , (54) 



where 5(i?, d) is the segment of a circle of radius d corresponding to a central angle 2$, 

5(1?, d) = d 2 h{$) , h{&) = i? - sin i? cos i? , and cos i? = V - . (55) 

CL 

Evaluating the integral one obtains 

q(L y ) = j- {*d 2 L y - U 3 ) (e~^ - 1) . (56) 
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FIG. 9: Interaction disks for the square well potential in Eq. (|13p for L y > d (top), and L y < d 
(bottom) in a channel with reflecting boundary condistions. 



For L y < d the overlap integral is given by (see the bottom panel of Fig. [U]) 

'* Li I, 



7id z L y - 2 



S($)dy 



1), 



(57) 



where, as in the previous case, S($) and y are defined by Eq. (155]) . The integral can be 
readily evaluated to yield 



1 



7T 



d 2 U 



2<i' { - + $ cos $ - sin i? + - sin 3 $ 
3 3 



-1) 



with 



cos-^r 



d 



(58) 



(59) 



The pressure components may be evaluated by using equations ( !56|58l) in the general ex- 
pressions ( 17151) . 

This analysis demonstrates that the pressure components exhibit a singularity at a chan- 
nel width L y = d. The nature of the singularity is obtained by expanding the overlap integral 
in 1 — Ly/d = e for small e > 0. In this limit one has $o — v2e and q(L y ) ~ — 2/15$q. Thus, 
for small e the singular part of the overlap integral Sq(L y ) satisfies 



Sq(Ly) 



r 5/2 



(60) 



20 
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FIG. 10: Channel- width dependence of the pressures for N = 60 particles, which interact with 
each other with the weak repulsive step potential of Eq. f)13[) . where u = 1. Reflecting boundary 
conditions are used as described in the main text. The particle diameter d = 1.5 in the reduced 
units applied, and the energy per particle, E/N, is unity. Keeping the particle density constant, 
N/V = 0.01, the temperature varies slightly with L y . 

As a result, P xx exhibits a 5/2 singularity at L y = d, and its third derivative with respect 
to L y diverges, as L y approaches d from below. On the other hand, P yy exhibits a stronger, 
3/2, singularity, as it is related to the derivative of the overlap integral with respect to L y . 

A comparison between theoretical and simulation results is provided in Fig. [TU] for a 
positive step potential. We use reduced units, in which the particle diameter d is 1.5, and 
for which the total energy per particle, E/N = Eq, is unity. In these units, we choose for 
the potential u — 1. In the simulation we studied N = 120 particles enclosed in a box 
with reflecting boundaries both parallel and perpendicular to the channel axis, such that 
the centers of the particles are confined to the volume V = L x L y (and the particle disks to 
the volume (L x + d)(L y +(£))■ Varying the channel width L y , the density is kept constant, 
N/V = 0.01. The singular width, d, is indicated by the vertical line. The agreement between 
the points from the simulation and the theoretical smooth lines corresponding to P xx , P yy 
and P = (P xx + Pyy)/ 2 is nearly perfect. 

To demonstrate the scaling directly, we plot in Fig. []J]the computer simulation results for 
the singular pressure contributions of P xx and P yy below the singular channel width L y = d. 
To do this, we note that the non-singular contribution to the overlap integral Eq. (!58l) for 
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FIG. 11: Channel- width scaling of the singular pressure contributions, which are obtained by 
subtracting the respective non-singular contributions from the computer simulation results for the 
pressures (points). The straight lines indicate the theoretically expected scaling. For details we 
refer to the main text. 

L y < d is given by Eq. (1561) . also evaluated at L y < d. The corresponding non-singular (NS) 
pressure contributions then follow from Eqs. (1201) and (|2"T1) with q(L y ) taken from Eq. (1561) . 
If this non-singular part is subtracted from the pressures determined by the simulations, a 
plot of 



PaaV 



- 1 



NS 



a e {x,y} 



as a function of the distance from the singularity, e = {d — L y )/d, reveals the expected 
scaling for the xx and yy pressure components. This is demonstrated by the straight lines 
in the log- log plot of Fig. [TTJ which are fully consistent with the expected scaling, 3/2 for 
P xx , and 5/2 for P yy . 

It is interesting to note that a system with a purely negative box potential, u < 0, is 
thermally unstable and tends to form clusters. The maximum entropy state consists of 
a single cluster of overlapping disks, floa ting in the gas of the remaining particles. As a 



consequence, the specific heat is negative jig, ll9|, |20j , and the temperature is increased due 



to the conversion of potential energy into kinetic. Such a property may only arise in the 
microcanonical ensemble and is familiar for gravitational systems. However, the attracting 
force need not be of long range 21J. A negative specific heat may even be observed for 



quantum-mechanical Coulomb systems 



20| and in experiments on nuclear fragmentation 
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FIG. 12: The interaction disks of a particle with a two-step potential in the case of reflecting 
boundaries. The inner disk is of diameter 2d and that of the outer disk is 2D. The boundaries of 
the channel are indicated by bold lines. The thin lines, are at a distance d/2 from the respective 
boundary, and they represent the limits, which the centers of the disks cannot cross. The vertical 
coordinate of the disk, y, and the angles (f> and $ are indicated. 



221 ] and atomic clusters 23j, |24| . In a preliminary study, we have observed the clustering 
for u — — 1, but we do not consider this case in more detail, because the times for reaching 
thermodynamic equilibrium are excessively long. 



B. Soft disks with a hard core: two-step potential 

In this section we consider particles which interact with a two-step potential 



u(r) 



oo r < d 

u d<r<D (61) 

r > D 



The soft potential u may be either attractive or repulsive. The disks interact with the walls 
of the channel only by the hard core interaction with diameter d (In the previous section this 
potential collapses into a 5 function). As before, we assume reflecting boundary conditions 
in the y direction as indicated in Figs. [12] and [T3], and periodic boundary conditions in 
the x direction. Thus, the volume accessible to the centers of the particles is given by 
V = (L y — d)L x . 

In calculating the overlap integral, one should distinguish between three regimes. 
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FIG. 13: The interaction disks of a particle with a two-step potential in the case of reflecting 
boundaries for L y < 2d, where the angles 4>o and $o are defined. In the intermediate regime 
2d < L y < D + d, the smaller disk does not intersect the boundary of the channel, and $o becomes 
0. 



For L y > D + d the integral may be expressed as (see Fig. (JT2 



q(L y ) = ir(D 2 - d 2 )^ - 1) - ;r 

rD 



d 2 



L y — d 



D z I h(4>)dy - d 2 / h{&)dy 
o Jo 



-/3u 



+ 



Ly — d 



d 2 / h(d)dy 



(62) 



where the function h is given in Eq. (tool) , and the angles and d are related to y via 
(see Fig. (TJJ) 

y = dcosfi = Dcos<f) . (63) 



Evaluating the integrals, one obtains 



q(L y 



4 D 3 



3 Ly d 



)(e~ /3u -1) - (nd 2 - - 



4 d 3 



-f3u 



3 L y — d 



For 2d < L y < D + d the overlap integral becomes 



(64) 



q(L y ) = n(D 2 - d 2 ){e~ pu - 1) - nd 
2 



Ly d 



+ 



D / h{(j))dy - d 2 / /i(i?)dy 

Z -2 



L y — d 



d l \ h{-&)dy 



o 



The integrals are readily evaluated to yield 



q(L y 



(ttD 



- (vrrf 2 



2 2iy ' g{<h)){<r*-i) 



Ly d 
4 d 3 , 
3 Ly — d' 



-flu 



-f3u 



1) 



(65) 



(66) 
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where 



and (see Fig. (fl3T 



/ N 2 ■ 1 • 3 

<7(a) = - + a cos a — sin a + - sin a , 



cos 0o 



(67) 
(68) 



For L y < 2d the overlap integral is 

q(L y ) = 7t(D 2 - d 2 ){e-? u 
2 



+ 



L y — d 



Ly — d 



7Td 2 



h((f))dy — d 2 



h($)dy 



-flu 



1) 



h(ti)dy . 



(69) 



It yields 



q{L y 



(nD 2 
(ird 2 - 



2D 3 



Ly — d 
2D 3 



#(0o))(e 



-/3u 



9(#o))e 



-/3u 



where 0o is given by Eq. 



and 



cos^n 



d 



L y — d 

d 



(70) 



(71) 



D + d and L y = 2d can 



(see Fig. (USD). 

The nature of the singularities of the pressure components at L y 
be analyzed as before. It is straightforward to show that at both points the P yy components 
exhibits a 3/2 singularity, while the P xx component exhibits a 5/2 singularity. 

Comparisons of these results with computer simulations for N = 60 particles of equal 
mass m are provided in Fig. [TH for the potential (16 ip with a positive step u = 1, and in 
Fig. [13 for the case of a negative step potential, u — — 1. All quantities are given in reduced 
units, for which the particle mass m, the hard core diameter, d, and the total energy per 
particle, E = E/N, are unity. The outer diameter is taken to be D = 1.5. The singular 
points at L y — d + D — 2.5 and L y = 2d = 2 are marked by the vertical lines. The density 
N/V = 0.01. For the computation of the theoretical pressures resulting in the smooth lines 
of Figs. dH and [13 the slight variation of the kinetic energy and, hence, of the temperature 
with the channel width was taken into account. The agreement between the theoretical 
expressions and the computer simulation results for the potential part of the pressures is 
very satisfactory. 
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FIG. 14: Channel width dependence of the potential-generated pressures for the interaction po- 
tential (|6ip with d = 1, D = 1.5 and u = 1. The dots are computer simulation results for iV = 60 
particles at a density ./V/V = 0.01. The smooth lines are obtained from the theoretical overlap 
integrals of Sec. IIVBI The singular channel widths are indicated by vertical lines. Reduced units 
are used, for which the particle mass m and the energy per particle, E/N, are unity. 
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FIG. 15: Channel width dependence of the potential-generated pressures for the interaction poten- 
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particles at a density N/V = 0.01. The smooth lines are obtained from the theoretical overlap 
integrals of Sec. IIVBI The singular channel widths are indicated by vertical lines. Reduced units 
are used, for which the particle mass m and the energy per particle, E/N, are unity. 
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V. SUMMARY 



In the paper we studied the pressure tensor of a system of disks moving in a narrow two 
dimensional channel, with either periodic or reflecting boundary conditions. We considered 
the low density regime using the Mayer cluster expansion, and tested the validity of the 
expansion using molecular dynamics studies. It is found that whenever the two-body inter- 
action potential between disks, u(r), exhibits a singularity at some distance ro, the pressure 
tensor exhibits a singularity as a function of the channel width, at one or more widths which 
are simply related to r . By studying several classes of interaction potentials, some rather 
general conclusions regarding the singularities of the pressure tensor can be reached. 

In the case of periodic boundary conditions, singularities take place at channel widths 
L y = 2ro/n with n = 1,2,... . For potentials which exhibit a discontinuity at ro, the 
transverse pressure, P yy , exhibits a 1/2 singularity while the longitudinal component, P xx , 
exhibits a weaker 3/2 singularity. For potentials which are continuous at r , and whose 
singular part vanishes as (ro — r) K , the transverse pressure exhibits a 2k + 1/2 singularity 
while the singularity of the longitudinal pressure is 2k + 3/2. Although these results have 
been demonstrated for specific interaction potentials u(r), they are rather general, as they 
are related only to the nature of the singularity of the potential. 

In the case of reflecting boundary conditions the pressure tensor exhibits a singularity at 
L y = r . The singularity is weaker than that of the case of periodic boundary conditions. 
Particularly, it was found that for a potential which is discontinuous at ro , the transverse 
component of the pressure exhibits a 3/2 singularity, while the longitudinal component 
exhibits a weaker 5/2 singularity. 



VI. APPENDIX 



There are a few minor misprints in some of the equations in the paper on hard disks 
[ijj] which are corrected below. These misprints do not affect any of the expressions for the 
pressure derived in that paper, or any of the numerical results. 
In particular, Equation (3) of Ref. [14( should read 



s 



S 
N 



In ( v — 



27 



and Equation (13) in [14| should become 

PxxV ^ f ds 



kT Lx ydL, 
P yy v _ L ( ds 



kT ' y \dL y/l 

Another misprint concerns the definition of v c in the expression for the virial in Eq. (8). v c 
is the velocity change of a particle % taking part in a binary collision c. 
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